% Display app res, phase & asso errors
% with impedance tensor & erros as input
% tested 9/11/02 -ok
% changed the error bar plotting 12/11/02
% known bug - if (Rhxy - RhxyErr/2) <0, log plot will propmt an error message
% latest date 13.10.2002

function[Rhxy,Rhyx,Phxy,Phyx] = cmpresph(tf1,ProcDef1,tf2,ProcDef2,cpp),
% cpp 1 rhoxy 2 rhoyx 3 rhoxy rhoyx phxy phyx
minT = 1/1000;
maxT = 120;

frq=ProcDef1.TLFreq1;

L = 1./frq > minT & 1./frq < maxT;

frq(~L) = [];
tf1(~L,:) = [];

file=ProcDef1.file1;
nfrq = length(frq);

for i = 1:nfrq,   
   Phxy(i) = atan2(imag(tf1(i,1)),real(tf1(i,1))) * 180 / 3.1415926535897;
   Rhxy(i) = 0.2 * (1/frq(i)) * (tf1(i,1) * conj(tf1(i,1)));
   Phyx(i) = (atan2(imag(tf1(i,2)),real(tf1(i,2))) * 180 / 3.1415926535897);
   Rhyx(i) = 0.2 * (1/frq(i)) * (tf1(i,2) * conj(tf1(i,2)));
end;

k=find(Phyx>-180&Phyx<-90);
Phyx(k)=Phyx(k)+180;
k=find(Phxy<-180&Phxy<-90);
Phxy(k)=Phxy(k)+180;

k=find(Phyx>-90&Phyx<-0);
Phyx(k)=Phyx(k)+90;
k=find(Phxy>-90&Phxy<-0);
Phxy(k)=Phxy(k)+90;


PhxyErr = (tf1(:,5)./abs(tf1(:,1)))*(180/pi); 
PhyxErr = (tf1(:,6)./abs(tf1(:,2)))*(180/pi);
RhxyErr = (tf1(:,5).*abs(tf1(:,1))./frq') * 0.4; % bit confusion. may be ok -dt 2.8.2
RhyxErr = (tf1(:,6).*abs(tf1(:,2))./frq') * 0.4;

%--------------------------------------------------------------
T = 1./frq;

mnx=ceil(log10(max(T)))+1;
mux=floor(log10(min(T)))-1;
mny=ceil(log10(max([Rhxy Rhyx])))+1;
muy=floor(log10(min([Rhyx Rhxy])))-1;

c1 = 'b.';
c2 = 'b.';
c3 = 'b';
c4 = 'b';


%subplot(2,1,1);

   plot(log10(T'),log10(Rhxy),c1);%,log10(T),log10(Rhyx),c2);
for i = 1:nfrq,
   tt = log10(Rhxy(i));
   kk = log10(Rhxy(i)+RhxyErr(i)/2);
   vv = 2*tt-kk;
   h=line(log10([T(i) T(i)]),[kk vv]);
   set(h,'color',c3);
end;

axis([mux mnx muy mny]);
grid;
ylabel('log ohm-m');
xlabel('log period (s)');
hold on;
% subplot(2,1,2);
% plot(log10(T),Phxy,c1);%,log10(T),Phyx,c2);
% 
% for i = 1:nfrq,
% 	h=line(log10([T(i) T(i)]),[Phxy(i)+PhxyErr(i)/2,Phxy(i)-PhxyErr(i)/2]);
%    set(h,'color',c3);
% end;

% axis([mux mnx 0 90]);
grid;

% xlabel('log period (s)');
% ylabel('degrees');
% title(file);
% hold on;

clear Rhxy Rhyx Phxy Phyx RhxyErr PhxyErr RhyxErr PhyxErr;

frq=ProcDef2.TLFreq1;

L = 1./frq > minT & 1./frq < maxT;

frq(~L) = [];
tf2(~L,:) = [];

file2=ProcDef2.file1;
nfrq = length(frq);

for i = 1:nfrq,   
   Phxy(i) = atan2(imag(tf2(i,1)),real(tf2(i,1))) * 180 / 3.1415926535897;
   Rhxy(i) = 0.2 * (1/frq(i)) * (tf2(i,1) * conj(tf2(i,1)));
   Phyx(i) = (atan2(imag(tf2(i,2)),real(tf2(i,2))) * 180 / 3.1415926535897);
   Rhyx(i) = 0.2 * (1/frq(i)) * (tf2(i,2) * conj(tf2(i,2)));
end;

k=find(Phyx>-180&Phyx<-90);
Phyx(k)=Phyx(k)+180;
k=find(Phxy<-180&Phxy<-90);
Phxy(k)=Phxy(k)+180;

k=find(Phyx>-90&Phyx<-0);
Phyx(k)=Phyx(k)+90;
k=find(Phxy>-90&Phxy<-0);
Phxy(k)=Phxy(k)+90;


PhxyErr = (tf2(:,5)./abs(tf2(:,1)))*(180/pi); 
PhyxErr = (tf2(:,6)./abs(tf2(:,2)))*(180/pi);
RhxyErr = (tf2(:,5).*abs(tf2(:,1))./frq') * 0.4; % bit confusion. may be ok -dt 2.8.2
RhyxErr = (tf2(:,6).*abs(tf2(:,2))./frq') * 0.4;

%--------------------------------------------------------------
T = 1./frq;

mnx=ceil(log10(max(T)))+1;
mux=floor(log10(min(T)))-1;
mny=ceil(log10(max([Rhxy Rhyx])))+1;
muy=floor(log10(min([Rhyx Rhxy])))-1;


%subplot(2,1,1);

if ishold == 1,
   c1 = 'r.';
   c2 = 'r.';
   c3 = 'red';
   c4 = 'red';
end;

   plot(log10(T'),log10(Rhxy),c1);%,log10(T),log10(Rhyx),c2);
for i = 1:nfrq,
   tt = log10(Rhxy(i));
   kk = log10(Rhxy(i)+RhxyErr(i)/2);
   vv = 2*tt-kk;
   h=line(log10([T(i) T(i)]),[kk vv]);
   set(h,'color',c3);
end;
axis([mux mnx muy mny]);
grid;
ylabel('log ohm-m');
xlabel('log period (s)');
%legend(file,'',file2);
hold on;
% subplot(2,1,2);
% plot(log10(T),Phxy,c1);%,log10(T),Phyx,c2);
% 
% for i = 1:nfrq,
% 	h=line(log10([T(i) T(i)]),[Phxy(i)+PhxyErr(i)/2,Phxy(i)-PhxyErr(i)/2]);
%    set(h,'color',c3);
% end;
% 
% axis([mux mnx 0 90]);
% grid;
% 
% xlabel('log period (s)');
% ylabel('degrees');
% %legend(file,'',file2);
% title(file);
% hold on;

